Introducción al Uso de R para Análisis Estadístico

Flujo de trabajo

Importar

Con R podemos acceder a información desde distintas fuentes:

  • Archivos de texto plano: txt, csv, tsv
  • Archivos MS excel, SPSS, SAS
  • Base de Datos: MySQL, SQLServer, PostgreSQL
  • APIs
library(tidyverse)

comunas <- read_csv("data/codigos_comunales.csv")
comunas
## # A tibble: 347 x 2
##    CODIGO COMUNA       
##     <int> <chr>        
##  1  15101 Arica        
##  2  15102 Camarones    
##  3  15201 Putre        
##  4  15202 General Lagos
##  5   1101 Iquique      
##  6   1402 Camiña       
##  7   1403 Colchane     
##  8   1404 Huara        
##  9   1405 Pica         
## 10   1401 Pozo Almonte 
## # ... with 337 more rows
library(haven)

casen <- read_sav("data/casen/Casen 2015.sav")
casen
## # A tibble: 266,968 x 776
##       folio     o id_vivienda region  provincia comuna zona  hogar tot_hog
##       <dbl> <dbl>       <dbl> <dbl+l> <dbl+lbl> <dbl+> <dbl> <dbl>   <dbl>
##  1  1.10e10     1  1101100101 1       11        1101   1         1       1
##  2  1.10e10     2  1101100101 1       11        1101   1         1       1
##  3  1.10e10     3  1101100101 1       11        1101   1         1       1
##  4  1.10e10     4  1101100101 1       11        1101   1         1       1
##  5  1.10e10     5  1101100101 1       11        1101   1         1       1
##  6  1.10e10     1  1101100102 1       11        1101   1         1       2
##  7  1.10e10     1  1101100103 1       11        1101   1         1       1
##  8  1.10e10     1  1101100104 1       11        1101   1         1       1
##  9  1.10e10     2  1101100104 1       11        1101   1         1       1
## 10  1.10e10     1  1101100105 1       11        1101   1         1       1
## # ... with 266,958 more rows, and 767 more variables: tot_par <dbl>,
## #   tot_nuc <dbl>, tot_per <dbl>, pres <dbl+lbl>, marca <chr>,
## #   pco1 <dbl+lbl>, sexo <dbl+lbl>, edad <dbl>, ecivil <dbl+lbl>,
## #   h5l <dbl>, h5h <dbl>, pareja <dbl+lbl>, nucleo <dbl+lbl>,
## #   pco2 <dbl+lbl>, e1 <dbl+lbl>, e2a <dbl+lbl>, e2b <dbl+lbl>,
## #   e2c <dbl+lbl>, e3 <dbl+lbl>, e4 <dbl+lbl>, e4_esp <chr>,
## #   e5a <dbl+lbl>, e5a_esp <chr>, e5b <dbl+lbl>, e6a <dbl+lbl>,
## #   e6b <dbl+lbl>, e6c_cod <dbl+lbl>, e6c <chr>, e6d <dbl+lbl>,
## #   e7nom <chr>, e7dir <chr>, e7com <chr>, e7rbd <dbl+lbl>, e7dv <dbl>,
## #   e7depen <dbl+lbl>, e7te <dbl+lbl>, e7rbd_sup <dbl+lbl>,
## #   e7sup_gratuidad <dbl+lbl>, e7com_cod <dbl+lbl>, e8 <dbl+lbl>,
## #   e9 <dbl+lbl>, e10a <dbl+lbl>, e10b <dbl+lbl>, e10c <dbl+lbl>,
## #   e10d <dbl+lbl>, e10e <dbl+lbl>, e11a <dbl+lbl>, e11bt1 <dbl+lbl>,
## #   e11bt1_esp <chr>, e11bt2 <dbl+lbl>, e11bt2_esp <chr>,
## #   e12pbu <dbl+lbl>, e12pbt <dbl+lbl>, e12pbd <dbl+lbl>,
## #   e12pbm <dbl+lbl>, e12bu <dbl+lbl>, e12bt <dbl+lbl>, e12bd <dbl+lbl>,
## #   e12bm <dbl+lbl>, e12bpc <dbl+lbl>, e12mu <dbl+lbl>, e12mt <dbl+lbl>,
## #   e12md <dbl+lbl>, e12mm <dbl+lbl>, e13a <dbl+lbl>, e13b <dbl+lbl>,
## #   e14 <dbl+lbl>, e15a <dbl+lbl>, e15b <dbl>, e16t1 <dbl+lbl>,
## #   e16t2 <dbl+lbl>, e0 <dbl+lbl>, o1 <dbl+lbl>, o2 <dbl+lbl>,
## #   o3 <dbl+lbl>, o4 <dbl+lbl>, o5 <dbl+lbl>, o6 <dbl+lbl>,
## #   o7r1 <dbl+lbl>, o7r2 <dbl+lbl>, o8 <dbl+lbl>, o9a <chr>, o9b <chr>,
## #   oficio1 <dbl+lbl>, oficio4 <chr+lbl>, o10 <dbl+lbl>, o11 <dbl+lbl>,
## #   o12 <dbl+lbl>, o13 <dbl+lbl>, o14 <dbl+lbl>, o15 <dbl+lbl>,
## #   o16 <dbl+lbl>, o17 <dbl+lbl>, o18 <dbl+lbl>, o19 <dbl+lbl>,
## #   o20 <dbl+lbl>, o21 <chr>, rama4_sub <chr+lbl>, rama1_sub <dbl+lbl>,
## #   o22 <chr>, ...
library(dplyr)
library(DBI)

con <- dbConnect(
  RMySQL::MySQL(),
  dbname = "censo2017",
  host = "142.93.20.188", 
  port = 3306,
  user = "test",
  password = "HFW9KYZBnEYr!"
)

dbListTables(con)
## [1] "hogar"    "idgeo"    "personas" "vivienda"
personas <- tbl(con,"personas")
personas
## # Source:   table<personas> [?? x 42]
## # Database: mysql 5.7.23-0ubuntu0.16.04.1 [test@142.93.20.188:/censo2017]
##    REGION PROVINCIA COMUNA    DC  AREA ZC_LOC ID_ZONA_LOC  NVIV NHOGAR
##     <int>     <int>  <int> <int> <int>  <int>       <int> <int>  <int>
##  1     15       152  15202     1     2      6       13225     1      1
##  2     15       152  15202     1     2      6       13225     3      1
##  3     15       152  15202     1     2      6       13225     3      1
##  4     15       152  15202     1     2      6       13225     3      1
##  5     15       152  15202     1     2      6       13225     3      1
##  6     15       152  15202     1     2      6       13225     9      1
##  7     15       152  15202     1     2      6       13225     9      1
##  8     15       152  15202     1     2      6       13225     9      1
##  9     15       152  15202     1     2      6       13225     9      1
## 10     15       152  15202     1     2      6       13225    10      1
## # ... with more rows, and 33 more variables: PERSONAN <int>, P07 <int>,
## #   P08 <int>, P09 <int>, P10 <int>, P10COMUNA <int>, P10PAIS <int>,
## #   P11 <int>, P11COMUNA <int>, P11PAIS <int>, P12 <int>, P12COMUNA <int>,
## #   P12PAIS <int>, P12A_LLEGADA <int>, P12A_TRAMO <int>, P13 <int>,
## #   P14 <int>, P15 <int>, P15A <int>, P16 <int>, P16A <int>,
## #   P16A_OTRO <int>, P17 <int>, P18 <chr>, P19 <int>, P20 <int>,
## #   P21M <int>, P21A <int>, P10PAIS_GRUPO <int>, P11PAIS_GRUPO <int>,
## #   P12PAIS_GRUPO <int>, ESCOLARIDAD <int>, P16A_GRUPO <int>

Transformar

casen_comuna <- casen %>% 
  mutate(comuna = as.numeric(comuna)) %>% 
  group_by(comuna) %>% 
  summarise(ingreso_promedio_mm = mean(y1, na.rm = TRUE)/1000)
casen_comuna
## # A tibble: 324 x 2
##    comuna ingreso_promedio_mm
##     <dbl>               <dbl>
##  1   1101                551.
##  2   1107                418.
##  3   1401                423.
##  4   1402                260.
##  5   1404                316.
##  6   1405                426.
##  7   2101                535.
##  8   2102                568.
##  9   2103                374.
## 10   2104                380.
## # ... with 314 more rows
personas_resumen <- personas %>% 
  group_by(region, comuna) %>% 
  summarise(personas = n(), escolaridad_promedio = mean(ESCOLARIDAD)) %>% 
  collect()

personas_resumen
## # A tibble: 346 x 4
## # Groups:   region [15]
##    region comuna personas escolaridad_promedio
##     <int>  <int>    <dbl>                <dbl>
##  1      1   1101   191468                 12.6
##  2      1   1107   108375                 12.2
##  3      1   1401    15711                 11.8
##  4      1   1402     1250                 15.2
##  5      1   1403     1728                 15.9
##  6      1   1404     2730                 13.1
##  7      1   1405     9296                 14.4
##  8      2   2101   361873                 13.3
##  9      2   2102    13467                 13.4
## 10      2   2103    10186                 13.9
## # ... with 336 more rows
data <- comunas %>%
  inner_join(personas_resumen, by = c("CODIGO" = "comuna")) %>% 
  inner_join(casen_comuna, by = c("CODIGO" = "comuna"))

data
## # A tibble: 324 x 6
##    CODIGO COMUNA        region personas escolaridad_prom~ ingreso_promedi~
##     <dbl> <chr>          <int>    <dbl>             <dbl>            <dbl>
##  1  15101 Arica             15   221364              12.9             402.
##  2  15102 Camarones         15     1255              11.4             314.
##  3  15201 Putre             15     2765              12.1             326.
##  4   1101 Iquique            1   191468              12.6             551.
##  5   1402 Camiña             1     1250              15.2             260.
##  6   1404 Huara              1     2730              13.1             316.
##  7   1405 Pica               1     9296              14.4             426.
##  8   1401 Pozo Almonte       1    15711              11.8             423.
##  9   1107 Alto Hospicio      1   108375              12.2             418.
## 10   2101 Antofagasta        2   361873              13.3             535.
## # ... with 314 more rows
data <- data %>% 
  mutate(
    region = factor(region),
    region2 = fct_lump(region, n = 8)
    )

data
## # A tibble: 324 x 7
##    CODIGO COMUNA region personas escolaridad_pro~ ingreso_promedi~ region2
##     <dbl> <chr>  <fct>     <dbl>            <dbl>            <dbl> <fct>  
##  1  15101 Arica  15       221364             12.9             402. Other  
##  2  15102 Camar~ 15         1255             11.4             314. Other  
##  3  15201 Putre  15         2765             12.1             326. Other  
##  4   1101 Iquiq~ 1        191468             12.6             551. Other  
##  5   1402 Camiña 1          1250             15.2             260. Other  
##  6   1404 Huara  1          2730             13.1             316. Other  
##  7   1405 Pica   1          9296             14.4             426. Other  
##  8   1401 Pozo ~ 1         15711             11.8             423. Other  
##  9   1107 Alto ~ 1        108375             12.2             418. Other  
## 10   2101 Antof~ 2        361873             13.3             535. Other  
## # ... with 314 more rows

Visualizar

  • Visualización: Principal caraterística
  • Paquete ggplot2. Muy poderoso.
  • Existen librerías interactivas

Graficar ingreso promedio y escolaridad promedio por comuna

p <- ggplot(data) +
  geom_point(aes(x = ingreso_promedio_mm, y = escolaridad_promedio,
                 label = COMUNA))

Agregamos más información

p <- ggplot(data) +
  geom_point(aes(x = ingreso_promedio_mm, y = escolaridad_promedio,
                 label = COMUNA, color = region, size= personas))

Detalles

p <- ggplot(data) +
  geom_point(aes(x = ingreso_promedio_mm, y = escolaridad_promedio,
                 label = COMUNA, color = region, size= personas),
             alpha = 0.75) +
  scale_color_viridis_d(option = "magma") +
  scale_x_continuous(trans = "log", labels = scales::comma,
                     breaks = seq(0, 1e3, by = 250))

Aquí parte la magia

p <- ggplot(data) +
  geom_point(aes(x = ingreso_promedio_mm, y = escolaridad_promedio,
                 label = COMUNA, color = region, size= personas),
             alpha = 0.75) +
  scale_color_viridis_d(option = "magma") +
  scale_x_continuous(trans = "log", labels = scales::comma,
                     breaks = seq(0, 1e3, by = 250)) +
  facet_wrap(~region2)

Aquí sigue la magia

p <- ggplot(data) +
  geom_point(aes(x = ingreso_promedio_mm, y = escolaridad_promedio,
                 label = COMUNA, color = region, size= personas),
             alpha = 0.75) +
  scale_color_viridis_d(option = "magma") +
  scale_x_continuous(trans = "log", labels = scales::comma,
                     breaks = seq(0, 1e3, by = 250)) +
  facet_wrap(~region2) +
  geom_smooth(aes(x = ingreso_promedio_mm, y = escolaridad_promedio),
              method = "lm", se = FALSE, color = "red", size = 1.2)

Importar + Transformar + Visualizar

library(sf)

dgeo <- read_sf("data/R13/Comuna.shp", layer = "Comuna")
dgeo
## Simple feature collection with 52 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -71.71523 ymin: -34.29093 xmax: -69.76999 ymax: -32.92194
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs
## # A tibble: 52 x 7
##    REGION PROVINCIA COMUNA DESC_REGIO               DESC_PROVI DESC_COMUN 
##    <chr>  <chr>     <chr>  <chr>                    <chr>      <chr>      
##  1 13     131       13114  REGIÓN METROPOLITANA DE~ SANTIAGO   LAS CONDES 
##  2 13     131       13115  REGIÓN METROPOLITANA DE~ SANTIAGO   LO BARNECH~
##  3 13     131       13132  REGIÓN METROPOLITANA DE~ SANTIAGO   VITACURA   
##  4 13     131       13107  REGIÓN METROPOLITANA DE~ SANTIAGO   HUECHURABA 
##  5 13     131       13124  REGIÓN METROPOLITANA DE~ SANTIAGO   PUDAHUEL   
##  6 13     131       13125  REGIÓN METROPOLITANA DE~ SANTIAGO   QUILICURA  
##  7 13     131       13102  REGIÓN METROPOLITANA DE~ SANTIAGO   CERRILLOS  
##  8 13     131       13116  REGIÓN METROPOLITANA DE~ SANTIAGO   LO ESPEJO  
##  9 13     131       13119  REGIÓN METROPOLITANA DE~ SANTIAGO   MAIPÚ      
## 10 13     131       13105  REGIÓN METROPOLITANA DE~ SANTIAGO   EL BOSQUE  
## # ... with 42 more rows, and 1 more variable: geometry <POLYGON [°]>

Importar + Transformar + Visualizar

ggplot() +
  geom_sf(data = dgeo) 

Importar + Transformar + Visualizar

library(classInt)
niveles <- c("bajo", "medio", "alto")

dgeo <- dgeo %>% 
  left_join(data, by = c("COMUNA" = "CODIGO")) %>% 
  mutate(
    escolaridad = classint(escolaridad_promedio, n = 3, style = "kmeans", labels = niveles), 
    ingreso = classint(ingreso_promedio_mm, n = 3, style = "kmeans", labels = niveles)
  )

glimpse(dgeo)
## Observations: 52
## Variables: 15
## $ REGION               <chr> "13", "13", "13", "13", "13", "13", "13",...
## $ PROVINCIA            <chr> "131", "131", "131", "131", "131", "131",...
## $ COMUNA               <dbl> 13114, 13115, 13132, 13107, 13124, 13125,...
## $ DESC_REGIO           <chr> "REGIÓN METROPOLITANA DE SANTIAGO", "REGI...
## $ DESC_PROVI           <chr> "SANTIAGO", "SANTIAGO", "SANTIAGO", "SANT...
## $ DESC_COMUN           <chr> "LAS CONDES", "LO BARNECHEA", "VITACURA",...
## $ COMUNA.y             <chr> "Las Condes", "Lo Barnechea", "Vitacura",...
## $ region               <fct> 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 1...
## $ personas             <dbl> 294838, 105833, 85384, 98671, 230293, 210...
## $ escolaridad_promedio <dbl> 15.3671, 12.9681, 14.8505, 12.9317, 12.38...
## $ ingreso_promedio_mm  <dbl> 1207.6004, 961.6909, 1388.5933, 366.3252,...
## $ region2              <fct> 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 1...
## $ escolaridad          <fct> alto, bajo, medio, bajo, bajo, bajo, bajo...
## $ ingreso              <fct> alto, alto, alto, bajo, medio, medio, baj...
## $ geometry             <POLYGON [m]> POLYGON ((362355.2 6307345,..., P...

Importar + Transformar + Visualizar

p2 <- ggplot() +
  geom_sf(data = select(dgeo, COMUNA, geometry),
          fill = "gray95", color = "gray80", size = 0.1) +
  geom_sf(data = dgeo, fill = "darkred", color = "gray80", size = 0.1) +
  scale_fill_viridis_d() +
  facet_grid(ingreso ~ escolaridad) +
  theme_minimal() 

Importar + Transformar + Visualizar

Modelar

Recordemos

Recordemos

data
## # A tibble: 324 x 7
##    CODIGO COMUNA region personas escolaridad_pro~ ingreso_promedi~ region2
##     <dbl> <chr>  <fct>     <dbl>            <dbl>            <dbl> <fct>  
##  1  15101 Arica  15       221364             12.9             402. Other  
##  2  15102 Camar~ 15         1255             11.4             314. Other  
##  3  15201 Putre  15         2765             12.1             326. Other  
##  4   1101 Iquiq~ 1        191468             12.6             551. Other  
##  5   1402 Camiña 1          1250             15.2             260. Other  
##  6   1404 Huara  1          2730             13.1             316. Other  
##  7   1405 Pica   1          9296             14.4             426. Other  
##  8   1401 Pozo ~ 1         15711             11.8             423. Other  
##  9   1107 Alto ~ 1        108375             12.2             418. Other  
## 10   2101 Antof~ 2        361873             13.3             535. Other  
## # ... with 314 more rows

Supongamos que quisieramos saber el efecto de la escolaridad por ingreso en comunas por cada region

datag <- data %>% 
  group_by(region2) %>% 
  nest()

datag
## # A tibble: 9 x 2
##   region2 data             
##   <fct>   <list>           
## 1 Other   <tibble [47 x 6]>
## 2 4       <tibble [15 x 6]>
## 3 5       <tibble [36 x 6]>
## 4 6       <tibble [33 x 6]>
## 5 7       <tibble [30 x 6]>
## 6 8       <tibble [54 x 6]>
## 7 9       <tibble [32 x 6]>
## 8 10      <tibble [25 x 6]>
## 9 13      <tibble [52 x 6]>
library(broom)

datag <- datag %>% 
  mutate(
    modelo = map(data, lm, formula = escolaridad_promedio ~ ingreso_promedio_mm),
    parametros = map(modelo, tidy)
  )

datag
## # A tibble: 9 x 4
##   region2 data              modelo   parametros          
##   <fct>   <list>            <list>   <list>              
## 1 Other   <tibble [47 x 6]> <S3: lm> <data.frame [2 x 5]>
## 2 4       <tibble [15 x 6]> <S3: lm> <data.frame [2 x 5]>
## 3 5       <tibble [36 x 6]> <S3: lm> <data.frame [2 x 5]>
## 4 6       <tibble [33 x 6]> <S3: lm> <data.frame [2 x 5]>
## 5 7       <tibble [30 x 6]> <S3: lm> <data.frame [2 x 5]>
## 6 8       <tibble [54 x 6]> <S3: lm> <data.frame [2 x 5]>
## 7 9       <tibble [32 x 6]> <S3: lm> <data.frame [2 x 5]>
## 8 10      <tibble [25 x 6]> <S3: lm> <data.frame [2 x 5]>
## 9 13      <tibble [52 x 6]> <S3: lm> <data.frame [2 x 5]>

Volver al mundo normal

dmods <- datag %>% 
  select(region2, parametros) %>% 
  unnest()

dmods
## # A tibble: 18 x 6
##    region2 term                 estimate std.error statistic  p.value
##    <fct>   <chr>                   <dbl>     <dbl>     <dbl>    <dbl>
##  1 Other   (Intercept)         11.5       0.798       14.4   1.85e-18
##  2 Other   ingreso_promedio_mm  0.00151   0.00192      0.786 4.36e- 1
##  3 4       (Intercept)          9.05      1.65         5.47  1.08e- 4
##  4 4       ingreso_promedio_mm  0.00599   0.00507      1.18  2.59e- 1
##  5 5       (Intercept)          9.93      0.870       11.4   3.53e-13
##  6 5       ingreso_promedio_mm  0.00698   0.00261      2.67  1.14e- 2
##  7 6       (Intercept)          8.11      1.00         8.07  4.08e- 9
##  8 6       ingreso_promedio_mm  0.00882   0.00331      2.66  1.22e- 2
##  9 7       (Intercept)          8.08      0.826        9.77  1.59e-10
## 10 7       ingreso_promedio_mm  0.00871   0.00302      2.88  7.46e- 3
## 11 8       (Intercept)          9.30      0.598       15.5   3.40e-21
## 12 8       ingreso_promedio_mm  0.00438   0.00191      2.30  2.58e- 2
## 13 9       (Intercept)          8.63      1.54         5.59  4.36e- 6
## 14 9       ingreso_promedio_mm  0.00859   0.00515      1.67  1.06e- 1
## 15 10      (Intercept)         10.6       1.45         7.31  1.93e- 7
## 16 10      ingreso_promedio_mm  0.000535  0.00468      0.114 9.10e- 1
## 17 13      (Intercept)         11.6       0.270       43.0   3.83e-41
## 18 13      ingreso_promedio_mm  0.00279   0.000514     5.44  1.61e- 6

Transformar

dmods <- dmods %>% 
  select(region2, term, estimate) %>% 
  spread(term, estimate) 

dmods
## # A tibble: 9 x 3
##   region2 `(Intercept)` ingreso_promedio_mm
##   <fct>           <dbl>               <dbl>
## 1 4                9.05            0.00599 
## 2 5                9.93            0.00698 
## 3 6                8.11            0.00882 
## 4 7                8.08            0.00871 
## 5 8                9.30            0.00438 
## 6 9                8.63            0.00859 
## 7 10              10.6             0.000535
## 8 13              11.6             0.00279 
## 9 Other           11.5             0.00151

Visualizar

library(ggrepel)

p3 <- ggplot(dmods, aes(ingreso_promedio_mm, `(Intercept)`)) +
  geom_point(size = 7, alpha = 0.6, color = "gray60") +
  geom_text_repel(aes(label = region2), force = 20) +
  scale_x_continuous(limits = c(0, NA))

Comunicar

Comunicar no es solo

  • Presentar un par de números
  • Mostrar un gráfico

Es

  • Atraer al receptor, y facilitarle la lectura
  • Contar la historía, dar razones

R poseee muchas otras formas de comunicar:

  • Gráficos interactivos (htmlwidgets)
  • Creación de documentos (words, pdf, html)
  • Aplicaciones/dashboard webs
# library(highcharter)
# hchart(data, "point", hcaes(ingreso_promedio_mm, escolaridad_promedio, 
#                             group = region2, size = personas)) %>% 
#   hc_add_theme(hc_theme_smpl())
library(plotly)
ggplotly(p, height = 600) %>% 
  config(displayModeBar = FALSE)
library(leaflet)

pal <- colorNumeric("viridis", NULL)

l <- leaflet(dgeo) %>%
  addTiles() %>%
  addPolygons(
    stroke = TRUE, color = "#DADADA", weight = 1,
    smoothFactor = 0.3, fillOpacity = 0.75,
    fillColor = ~pal(escolaridad_promedio),
    label = ~paste0(DESC_COMUN, ": ", escolaridad_promedio)
    ) %>%
  addLegend(pal = pal, values = ~escolaridad_promedio, opacity = 1.0)
l

Gracias

Información Relevante: